home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
newmat03.lha
/
newmat03
/
newmat4.cxx
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
12KB
|
404 lines
//$$ newmat4.cxx Constructors, ReDimension, basic utilities
// Copyright (C) 1991: R B Davies and DSIR
#include "include.hxx"
#include "newmat.hxx"
#include "newmatrc.hxx"
//#define REPORT { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
#define REPORT {}
//#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
// REPORT1 constructors only - doesn't work in turbo and Borland C++
#define REPORT1 {}
//#define MONITOR(what,storage,store) \
// { cout << what << " " << storage << " at " << (long)store << "\n"; }
#define MONITOR(what,store,storage) {}
/*************************** general utilities *************************/
static int tristore(int n) // els in triangular matrix
{ return (n*(n+1))/2; }
/****************************** constructors ***************************/
GeneralMatrix::GeneralMatrix()
{ store=0; storage=0; nrows=0; ncols=0; tag=-1; }
GeneralMatrix::GeneralMatrix(int s)
{
REPORT1
storage=s; tag=-1;
store = new real [storage]; MatrixErrorNoSpace(store);
MONITOR("Make (GenMatrix)",storage,store)
}
Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
{ REPORT1 nrows=m; ncols=n; }
SymmetricMatrix::SymmetricMatrix(int n) : GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }
UpperTriangularMatrix::UpperTriangularMatrix(int n)
: GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }
LowerTriangularMatrix::LowerTriangularMatrix(int n)
: GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }
DiagonalMatrix::DiagonalMatrix(int m) : GeneralMatrix(m)
{ REPORT1 nrows=m; ncols=m; }
Matrix::Matrix(BaseMatrix& M)
{
REPORT1 CheckConversion(M);
GeneralMatrix* gmx=M.Evaluate(MatrixType::Rect); GetMatrix(gmx);
}
RowVector::RowVector(BaseMatrix& M) : Matrix(M)
{
if (nrows!=1) MatrixError("Illegal conversion to row vector"); }
ColumnVector::ColumnVector(BaseMatrix& M) : Matrix(M)
{
if (ncols!=1) MatrixError("Illegal conversion to column vector"); }
SymmetricMatrix::SymmetricMatrix(BaseMatrix& M)
{
REPORT1 CheckConversion(M);
GeneralMatrix* gmx=M.Evaluate(MatrixType::Sym); GetMatrix(gmx);
}
UpperTriangularMatrix::UpperTriangularMatrix(BaseMatrix& M)
{
REPORT1 CheckConversion(M);
GeneralMatrix* gmx=M.Evaluate(MatrixType::UT); GetMatrix(gmx);
}
LowerTriangularMatrix::LowerTriangularMatrix(BaseMatrix& M)
{
REPORT1 CheckConversion(M);
GeneralMatrix* gmx=M.Evaluate(MatrixType::LT); GetMatrix(gmx);
}
DiagonalMatrix::DiagonalMatrix(BaseMatrix& M)
{
REPORT1 CheckConversion(M);
GeneralMatrix* gmx=M.Evaluate(MatrixType::Diag); GetMatrix(gmx);
}
GeneralMatrix::~GeneralMatrix()
{
if (store)
{
MONITOR("Free (GenMatrix)",storage,store)
delete [storage] store;
}
}
CroutMatrix::CroutMatrix(BaseMatrix& m)
{
REPORT1
GeneralMatrix* gm = m.Evaluate(MatrixType::Rect); GetMatrix(gm);
if (nrows!=ncols) MatrixError("Matrix not square");
d=TRUE; indx=new int [nrows]; MatrixErrorNoSpace(indx);
ludcmp();
}
ReturnMatrix::ReturnMatrix(GeneralMatrix& gmx)
{
REPORT1
gm = gmx.Type().New(); MatrixErrorNoSpace(gm);
gm->GetMatrix(&gmx); gm->ReleaseAndDelete();
}
/**************************** ReDimension matrices ***************************/
void GeneralMatrix::ReDimension(int nr, int nc, int s)
{
REPORT
if (store)
{
MONITOR("Free (ReDimensi)",storage,store)
delete [storage] store;
}
storage=s; nrows=nr; ncols=nc; tag=-1;
store = new real [storage]; MatrixErrorNoSpace(store);
MONITOR("Make (ReDimensi)",storage,store)
}
void Matrix::ReDimension(int nr, int nc)
{ REPORT GeneralMatrix::ReDimension(nr,nc,nr*nc); }
void SymmetricMatrix::ReDimension(int nr)
{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
void UpperTriangularMatrix::ReDimension(int nr)
{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
void LowerTriangularMatrix::ReDimension(int nr)
{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
void DiagonalMatrix::ReDimension(int nr)
{ REPORT GeneralMatrix::ReDimension(nr,nr,nr); }
void RowVector::ReDimension(int nc)
{ REPORT GeneralMatrix::ReDimension(1,nc,nc); }
void ColumnVector::ReDimension(int nr)
{ REPORT GeneralMatrix::ReDimension(nr,1,nr); }
/********************* manipulate types, storage **************************/
int GeneralMatrix::search(const GeneralMatrix* s) const
{ REPORT return (s==this) ? 1 : 0; }
int AddedMatrix::search(const GeneralMatrix* s) const
{ REPORT return bm1->search(s) + bm2->search(s); }
int ShiftedMatrix::search(const GeneralMatrix* s) const
{ REPORT return bm->search(s); }
int NegatedMatrix::search(const GeneralMatrix* s) const
{ REPORT return bm->search(s); }
int ConstMatrix::search(const GeneralMatrix* s) const
{ REPORT return (s==cgm) ? 1 : 0; }
int ReturnMatrix::search(const GeneralMatrix* s) const
{ REPORT return (s==gm) ? 1 : 0; }
MatrixType Matrix::Type() const { return MatrixType::Rect; }
MatrixType SymmetricMatrix::Type() const { return MatrixType::Sym; }
MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
MatrixType DiagonalMatrix::Type() const { return MatrixType::Diag; }
MatrixType RowVector::Type() const { return MatrixType::RowV; }
MatrixType ColumnVector::Type() const { return MatrixType::ColV; }
MatrixType CroutMatrix::Type() const { return MatrixType::Crout; }
MatrixType RowedMatrix::Type() const { return MatrixType::RowV; }
MatrixType ColedMatrix::Type() const { return MatrixType::ColV; }
MatrixType DiagedMatrix::Type() const { return MatrixType::Diag; }
MatrixType MatedMatrix::Type() const { return MatrixType::Rect; }
MatrixType GetSubMatrix::Type() const { return mt; }
MatrixType AddedMatrix::Type() const
{ REPORT return bm1->Type() + bm2->Type(); }
MatrixType MultipliedMatrix::Type() const
{ REPORT return bm1->Type() * bm2->Type(); }
MatrixType SolvedMatrix::Type() const
{ REPORT return bm1->Type().i() * bm2->Type(); }
MatrixType SubtractedMatrix::Type() const
{ REPORT return bm1->Type() - bm2->Type(); }
MatrixType ShiftedMatrix::Type() const
{ REPORT MatrixType mteqel(MatrixType::EqEl); return bm->Type() + mteqel; }
MatrixType ScaledMatrix::Type() const { REPORT return bm->Type(); }
MatrixType TransposedMatrix::Type() const { REPORT return bm->Type().t(); }
MatrixType NegatedMatrix::Type() const { REPORT return bm->Type(); }
MatrixType InvertedMatrix::Type() const { REPORT return bm->Type().i(); }
MatrixType ConstMatrix::Type() const { REPORT return cgm->Type(); }
MatrixType ReturnMatrix::Type() const { REPORT return gm->Type(); }
int GeneralMatrix::NrowsV() const { return nrows; }
int RowedMatrix::NrowsV() const { return 1; }
int MatedMatrix::NrowsV() const { return nr; }
int GetSubMatrix::NrowsV() const { return row_number; }
int AddedMatrix::NrowsV() const { return bm1->NrowsV(); }
int ShiftedMatrix::NrowsV() const { return bm->NrowsV(); }
int TransposedMatrix::NrowsV() const { return bm->NcolsV(); }
int NegatedMatrix::NrowsV() const { return bm->NrowsV(); }
int ColedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
int DiagedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
int ConstMatrix::NrowsV() const { return cgm->Nrows(); }
int ReturnMatrix::NrowsV() const { return gm->Nrows(); }
int GeneralMatrix::NcolsV() const { return ncols; }
int ColedMatrix::NcolsV() const { return 1; }
int MatedMatrix::NcolsV() const { return nc; }
int GetSubMatrix::NcolsV() const { return col_number; }
int AddedMatrix::NcolsV() const { return bm2->NcolsV(); }
int ShiftedMatrix::NcolsV() const { return bm->NcolsV(); }
int TransposedMatrix::NcolsV() const { return bm->NrowsV(); }
int NegatedMatrix::NcolsV() const { return bm->NcolsV(); }
int RowedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int DiagedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int ConstMatrix::NcolsV() const { return cgm->Ncols(); }
int ReturnMatrix::NcolsV() const { return gm->Ncols(); }
// Rules regarding tDelete, reuse, GetStore
// All matrices processed during expression evaluation must be subject
// to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
// If reuse returns TRUE the matrix must be reused.
// GetMatrix(gm) always calls gm->GetStore()
// gm->Evaluate obeys rules
// bm->Evaluate obeys rules for matrices in bm structure
void GeneralMatrix::tDelete()
{
if (tag<0)
{
if (tag<-1) { REPORT store=0; delete this; return; } // borrowed
else { REPORT return; }
}
if (tag==1)
{
REPORT MONITOR("Free (tDelete)",storage,store)
if (store) delete [storage] store; store=0; tag=-1; return;
}
if (tag==0) { REPORT delete this; return; }
REPORT tag--; return;
}
BOOL GeneralMatrix::reuse()
{
if (tag<-1)
{
REPORT
real* s = new real [storage]; MatrixErrorNoSpace(s);
MONITOR("Make (reuse)",storage,s)
real* s1=store; real* s2=s; int i=storage;
while(i--) *s2++ = *s1++;
store=s; tag=0; return TRUE;
}
if (tag<0) { REPORT return FALSE; }
if (tag<=1) { REPORT return TRUE; }
REPORT tag--; return FALSE;
}
real* GeneralMatrix::GetStore()
{
if (tag<0 || tag>1)
{
real* s = new real [storage]; MatrixErrorNoSpace(s);
MONITOR("Make (GetStore)",storage,s)
real* s1=store; real* s2=s; int i=storage;
while(i--) *s2++ = *s1++;
if (tag>1) { REPORT tag--; }
else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
else { REPORT }
return s;
}
real* s=store; store=0;
if (tag==0) { REPORT delete this; }
else { REPORT tag=-1; }
return s;
}
#ifndef __ZTC__
void GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
{
REPORT tag=-1;
nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;
store = new real [storage]; MatrixErrorNoSpace(store);
MONITOR("Make (GetMatrix)",storage,store)
real* s1=gmx->store; real* s2=store; int i=storage;
while(i--) *s2++ = *s1++;
}
#endif
void GeneralMatrix::GetMatrix(GeneralMatrix* gmx)
{
REPORT tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
storage=gmx->storage; store=gmx->GetStore();
}
GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
// Copy storage of *this to storage of *gmx. Then convert to type mt.
// If mt == 0 just let *gm point to storage of *this if tag<0.
{
if ((int)mt==0)
{
if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
}
else if (mt!=gmx->Type())
{
REPORT gmx->tag = -2; gmx->store = store;
gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
}
else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
return gmx;
}
void GeneralMatrix::Eq(BaseMatrix& X, MatrixType mt)
// Count number of references to this in X.
// If zero delete storage in X;
// otherwise tag X to show when storage can be deleted
// evaluate X and copy to current object
{
int counter=X.search(this);
if (counter==0)
{
REPORT MONITOR("Free (operator=)",storage,store)
if (store) { REPORT delete [storage] store; storage=0; }
}
else { REPORT Release(counter); }
GeneralMatrix* gmx = X.Evaluate(mt);
if (gmx!=this) { REPORT GetMatrix(gmx); }
else { REPORT }
Protect();
}
void GeneralMatrix::Inject(const GeneralMatrix& X)
// copy stored values of X; otherwise leave els of *this unchanged
{
REPORT
CheckConversion(X);
if (nrows != X.nrows || ncols != X.ncols)
MatrixError("Inject: dimensions don't match");
MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
int i=nrows;
while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
}
void GeneralMatrix::CheckConversion(const BaseMatrix& M)
{ if (!(this->Type() >= M.Type())) MatrixError("Illegal Conversion"); }
/************************* nricMatrix routines *****************************/
void nricMatrix::MakeRowPointer()
{
row_pointer = new real* [nrows]; MatrixErrorNoSpace(row_pointer);
real* s = Store() - 1; int i = nrows; real** rp = row_pointer;
while (i--) { *rp++ = s; s+=ncols; }
}
void nricMatrix::DeleteRowPointer()
{ if (nrows) delete [nrows] row_pointer; }
void GeneralMatrix::CheckStore() const
{ if (!store) MatrixError("NRIC accessing matrix with dimensions not set"); }